# Replication for:
# Jesus was a Refugee: Religious Values Framing Can Increase Support for Refugees Among White Evangelical Republicans

# NOTE: Leave all of the file where they are, and this should run through the RProj included.
# If for some reason the data isn't loading properly, you may need to change your file path.
# setwd(here::here())
# or type manually to wherever the folder lives on your local device:
# setwd("")

source("00. Read Data.R", print.eval = F, echo = F)


# ---------------------------------------------------------------------------------------------------------------
# Table 1. Experimental Design and Treatment
# ---------------------------------------------------------------------------------------------------------------
df %>%
  tabyl(messages) %>%
  adorn_pct_formatting()

# ---------------------------------------------------------------------------------------------------------------

# Alpha for resettlement measures
x <- cbind(df$BCV2_recode, df$BCV3_recode)
psych::alpha(as.data.frame(x))

# Alpha for benefits measures
x <- cbind(df$BR1_recode, df$BR2_recode)
psych::alpha(as.data.frame(x))

# Alpha Evangelical identity
x <- cbind(df$R4_1_recode, df$R4_2_recode, df$R4_3_recode, df$R4_4_recode, df$R4_5_recode, df$R4_6_recode, df$R4_7_recode, df$R4_8_recode, df$R4_9_recode, df$R4_10_recode)
psych::alpha(as.data.frame(x))

# Alpha Republican identity
x <- cbind(df$e4_1_recode, df$e4_2_recode, df$e4_3_recode, df$e4_4_recode, df$e4_5_recode, df$e4_6_recode, df$e4_7_recode, df$e4_8_recode, df$e4_9_recode, df$e4_10_recode)
psych::alpha(as.data.frame(x),warnings = F)

# ---------------------------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------------------------
# Figure 1. Distribution of evangelical and Republican identities
# ---------------------------------------------------------------------------------------------------------------

df %>%
  dplyr::select(Republican_Index,weight_bornagain) %>%
  mutate(Var = 'Republican Index') %>%
  rename(Distribution = Republican_Index) -> RepsVar

df %>%
  dplyr::select(Evangelical_Index,weight_bornagain) %>%
  mutate(Var = 'Evangelical Index') %>%
  rename(Distribution = Evangelical_Index) -> EvangVar

MASTER <- rbind(RepsVar,EvangVar)

ggplot(data=MASTER, aes(x=Distribution, group=Var, fill=Var), weight = weight_bornagain) +
  geom_density(adjust=1.5, alpha=.3) +
  labs(title="Distibution of the Identity Strength Indices", 
       x="0 = Weak Identity to 1 = Strong Identity", y = "Density", caption = "", subtitle = "", fill = "Identity:") + 
  theme(legend.position = "none", axis.text.x = element_text(size=12, angle=35)) +
  scale_fill_grey(start = 0, end = 1) + 
  scale_color_grey() +
  steph_theme

# ---------------------------------------------------------------------------------------------------------------
# Table 2. Religious values message treatment effects 
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages, data = df, weights = weight_bornagain)
mod2 <- lm(FuturePolicy_Scale ~ messages, data = df, weights = weight_bornagain)
mod3 <- lm(BR1_recode ~ messages, data = df, weights = weight_bornagain)
mod4 <- lm(BR2_recode ~ messages, data = df, weights = weight_bornagain)
export_summs(mod1,mod2,mod3,mod4,
             error_format = "({std.error})",
             error_pos = c("below"),
             ci_level = 0.90,
             stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01),
             model.names = c('Feeling Therm.','Resettlement', 'Stay in School','Get Benefits'),
             coefs = NULL)

# ---------------------------------------------------------------------------------------------------------------
# Figure 2. Treatment effects (relative to control) by evangelical identity strength
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod3c <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod2c <- lm(BR1_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod4c <- lm(BR2_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
scaleFUN <- function(x) sprintf("%.3f", x)

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("RV+SC Interaction: Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("RV+SC Interaction: Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("RV+SC Interaction: Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("RV+SC Interaction: Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("RV Interaction: Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("RV Interaction: Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))


plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("RV Interaction: Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "Marginal effects of the Religious values message relative to the control group by the respondent's evangelical identification. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith evangelical identification. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("RV Interaction: Unemployment"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4)

# ---------------------------------------------------------------------------------------------------------------
# Figure 3. Treatment effects (relative to control) by Republican identity strength
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod3c <- lm(FuturePolicy_Scale ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod2c <- lm(BR1_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod4c <- lm(BR2_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("RV+SC Interaction: Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))


plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("RV+SC Interaction: Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("RV+SC Interaction: Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("RV+SC Interaction: Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("RV Interaction: Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))


plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("RV Interaction: Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))


plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("RV Interaction: Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "Marginal effects of the Religious values message relative to the control group by the respondent's Republican identification. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith Republican identification. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("RV Interaction: Unemployment"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4)


# ---------------------------------------------------------------------------------------------------------------
# Table 3. Religious values message treatment effects among non-evangelicals 
# ---------------------------------------------------------------------------------------------------------------
CCES %>%
  mutate(pew_bornagain_r = as_factor(pew_bornagain)) %>%
  filter(pew_bornagain_r == "No") -> CCES_nonevan

mod1 <- lm(FT1_1_recode ~ messages, data = CCES_nonevan, weights = teamweight)
mod2 <- lm(BR1_recode ~ messages, data = CCES_nonevan, weights = teamweight)
mod3 <- lm(BR2_recode ~ messages, data = CCES_nonevan, weights = teamweight)
mod6 <- lm(FuturePolicy_Scale ~ messages, data = CCES_nonevan, weights = teamweight)

jtools::export_summs(mod1,mod6,mod2,mod3,
             error_pos = c("below"),
             ci_level = 0.90,
             stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01),
             model.names = c('Feeling Therm.', 'Resettlement', 'Schooling','Unemployment'),
             coefs = NULL)





